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Abstract 

We discuss, in a pedagogical way, how to solve for relativistic wave functions from the radial 
Dirac equations. After an brief introduction, in Section II we solve the equations for a linear 
Lorentz scalar potential, V s (r), that provides for confinement of a quark. The case of massless 
u and d quarks is treated first, as these are necessarily quite relativistic. We use an iterative 
procedure to find the eigenenergies and the upper and lower component wave functions for the 
ground state and then, later, some excited states. Solutions for the massive quarks (s, c, and b) 
are also presented. In Section III we solve for the case of a Coulomb potential, which is a time- 
like component of a Lorentz vector potential, V v (r). We re-derive, numerically, the (analytically 
well-known) relativistic hydrogen atom eigenenergies and wave functions, and later extend that to 
the cases of heavier one-electron atoms and muonic atoms. Finally, Section IV finds solutions for 
a combination of the V s and V v potentials. We treat two cases. The first is one in which V s is the 
linear potential used in Sec. II and V v is Coulombic, as in Sec. III. The other is when both V s 
and V v are linearly confining, and we establish when these potentials give a vanishing spin-orbit 
interaction (as has been shown to be the case in quark models of the hadronic spectrum). 
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I. INTRODUCTION 

Dirac formulated relativistic quantum mechanics in the late 1920's.[l| Since field theory 
had not yet been developed, the relativistic aspect led to a number of confusions related 
to currents since only total charge is conserved. With the recognition of the appearance 
of antiparticles, and the development of field theory, analogous to the transition from the 
Canonical to the Grand Canonical ensemble in statistical mechanics, these issues were re- 
solved. Nonetheless, the solution of the equation itself, ignoring these deeper aspects, has 
proven valuable even in more modern contexts. 

For example, the MIT bag model of quark confinement, and the theories that evolved 
from it (chiral bag, cloudy bag, etc.) 3»] depend on solving the Dirac equation as a wave 
function for a particle in an effective potential. Through boundary conditions and other 
approximations, the Dirac equation has even been employed in this way in the study of 
nuclear structure. 0, 0] 

In particle physics, the quark model has been employed, with great success, to describe 
3aryon and meson states and their structure. Usually this is done in a non-relativistic model, 
6] where the potentials follow the patterns expected by a Foldy-Wouthuysen reduction [3] 
of a simpler, theoretically motivated potential in the Dirac equation. In our case, however, 
we wanted to see how these structures and solutions appear without the non-relativistic 
reduction approximations, by solving the Dirac equation itself for the simpler potential with 
fewer adjustable parameters. 

Going back somewhat in time, very soon after Dirac's initial formulation, Darwin js] 
found an analytic solution (!) for the radial wave functions for hydrogen-like atoms in terms 
of confluent hypergeometric functions. Since then, solutions of the Dirac equation have been 
important in atomic physics, even to recent times, jol, [lo| For example, more complicated 
central potentials than Za/r, the fourth component of a Lorentz four- vector, usually require 
finding a numerical solution. Another example is a muonic atom, in which a muon is in an 

n 

s-state about a heavy nucleus having a realistic charge distribution. 

Finding such numerical solutions involves solving coupled ordinary differential equations 
for the upper and lower components of the Dirac wave function, Q, Q 
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where I' = 2j — I and the y l j m are the spin-orbital angular momentum wave functions, 

yjm( e > 4>)=Y,< J m \ l h m ~ m s m s> Y l m{d, 4>) Xm a , (2) 

m s 

where < jm \ I |, m — m s m s > is a Clebsch-Gordan coefficient and x is a Pauli spinor. 

We consider here the case of bound-state wave functions for a central potential with 
energy eigenvalue E. For example, for hydrogen-like atoms the potential is V v (r) = —Za/r. 
With appropriate boundary conditions, the radial wave functions g(r) and f(r), which can be 
taken as real functions, are solutions of the following coupled first-order ordinary differential 
equations (ODEs), Q, Q 

9'{r) + ~0(r) - (E - V v (r) + m) f(r) = , (3) 
f'(r) - -g(r) + (E - K(r) - m) ^(r) = . (4) 

Here the integer k is determined by the angular momentum quantum numbers according to 

k = -{1 + 1) if j = l + \ , 

k = I if j = l~\ ■ (5) 

The reason for the subscript on the potential V v in those equations is to indicate that it 
is the fourth component of a Lorentz four-vector, such as the Coulomb potential. However, 
our motivation is not so much in atomic physics applications as it is for treating mesons as 
qq states in a relativistic quark model. The non-relativistic quark model often assumes the 
confining potential for quarks to be the so-called "Cornell potential," 6] 

V(r) = -— + K 2 r , (6) 
r 

where a s = g s 2 /4:7i, with g s being the (running) quark-gluon coupling constant, and k 2 (or 
a) is the string tension. It is the linear term in V(r) that confines the quarks, corresponding 
to the effective potential in the relativistic Bag Model. Q 

For a relativistic quark model, the two pieces of that non-relativistic potential have dif- 
ferent Lorentz transformation properties. The color Coulomb potential, a s /r, is the fourth 
component of a Lorentz vector, while the confining linear potential transforms as a Lorentz 
scalar, which we will write as V s (r). Thus, in solving the relativistic radial Dirac equations, 
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the two terms in the non-relativistic potential of Eq. fl6]) enter the coupled ODE's differently. 
The equations to be solved are 

g'(r) + *g(r) - (E - V v (r) + V s (r) + m) f(r) = , (7) 

f'(r) - -f(r) + (E - K(r) - V s (r) - m) g(r) = . (8) 
r 

Note that in these equations the sign for V v is opposite to that of the energy eigenvalue E, 
the fourth component of the momentum four-vector, while that for V s matches that for the 
fermion mass m, also a Lorentz scalar. 

So, we were motivated to study solutions of Eqs. ([7]) and (jHJ). After submitting an earlier 
version of this paper to the archives, we learned of related studies of these equations by 
Paris 1^] and by Soares de Castro and Franklin [ljj], who discussed analytic solutions. We, 
however, have adopted a numerical approach and, in this paper, we discuss how to go about 
solving these coupled ODE's numerically. The presentation follows our numerical journey 
more or less historically. Our hope in writing this paper is that, by reading it, students 
might avoid some of the pitfalls we encountered and then overcame. 

It turns out that solving for a linear scalar potential 0, [l^ is quite a bit easier (numer- 
ically) than the relativistic hydrogen atom, so we treat that first in Sec. II. We then turn, 
in Sec. Ill, to the problems we had with in solving the radial equations for a Coulomb-like 
potential for hydrogen-like atoms and how we dealt with them. This is followed, in Sec. IV, 
by a discussion of the mixed problem having both a scalar and a vector potential, as in Eqs. 
© and ©. 

II. LINEAR SCALAR POTENTIAL 



For the case of massless quarks (an approximation we have made 13 1 for the u and d 

nn 

quarks), we want to solve the equations for a confining linear scalar potential, |4J, [16j 

V„(r) = K 2 (r-r ) , (9) 

where (with h~ — c — 1), k has dimensions of fm _1 . The negative offset, —K 2 r , effectively 
gives these quarks a constituent mass. [4] It also provides a rough representation of the 
effect of how the short-distance color Coulomb interaction between quarks leads to quark 
confinement, albeit without the correct Lorentz representation properties. 



The case of massless quarks exhibits the role of relativity in a maximal way, meaning that 
the lower component f(x) is comparable in size to the upper component g(x). We will, later, 
choose an appropriate mass m for the massive quark flavors, s, c, and b, and solve for those 
wave functions. In these cases, as the quark mass m increases, f(x) becomes smaller relative 
to g(x), reflecting the transition to a non-relativistic limit where f(x) vanishes. Also, as will 
be seen, with increasing mass m, the heavier quark radial wave functions are progressively 
narrower than in the massless case. Thus we will, in Sec. IIVt also include the effects of a color 
Coulomb-like vector potential V v coming from one-gluon exchange. This should prevent a 
too-rapid fall-off of the quark-quark interaction energies as the quark masses increase 
but that is a matter for a separate paper. 

It simplifies the equations if we convert the ODEs to dimensionless form, dividing through 
by k and defining a dimensionless distance x — nr: 

g'(x) + ^g(x) - {E + V s (x) + m) f(x) = , (10) 
f'(x) - + (E — V s (x) - m) g(x) = , (11) 

where E = E/k, V s (x) = x — xq, and fh = mj 'k are now also dimensionless. 

Numerical integrations of these equations can be done by one or another of the Runge- 
Kutta methods. IH] To integrate these two first-order ODEs one needs to specify two bound- 
ary conditions (BCs) at the place where one begins the integration. However, as the radial 
wave functions eventually need to satisfy a normalization condition, 

POO POO 

/ x 2 dx tyl(x) + ^ b {x)} = / dx [g\x) + f(x)] = 1 , (12) 
Jo Jo 

this will determine one of the BCs, for example, the scale of the wave functions at asymp- 
totically large distances. 



A. Shooting Outwards 

In our first attempt to solve these equations we thought to integrate them outward from 
the origin with the hope of adjusting the energy eigenvalue E to assure that both g(x) and 
f(x) fall off to zero asymptotically. The first problem in doing this is the singularity in the 
equations at x = 0. This singularity, however, is easily avoided by starting instead at a 
small value of x = e away from the origin. 
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For starting values (BCs) we can take advantage of the expected power-law behaviors 
of ipa{x) and t/>&(se) near the origin for their known orbital angular momenta I and I'. For 
example, if j = I + ~, for which k = —(I + 1) and the lower component's I' — I + 1, one can 
choose 

g(e) = a e l+1 , f(e)=b e l+2 , (13) 

as motivated by the non-relativistic limit, where ip a — > \1/nr ~ r ' at small r. Substituting 
this /(e) into the second ODE allows us to determine bo: 

/(e) = -ao (E - V s (e) - m) e l+2 /(2l + 3) . (14) 

The parameter ao at this point is arbitrary for now and will be fixed later by the normal- 
ization condition, Eq. ( Fl2l) . 

The problem with simply shooting outwards is what to choose for the energy E. The 
ODE solver for almost all initial guesses for E will have the well-known problem that the 
calculated g(r) and /(r) very soon blow up exponentially, either positively or negatively. It 
might be possible to iteratively refine the guess for E to eventually find decaying solution, 
but we soon decided to try a different method. 



B. Shooting Inwards 

A better way of assuring decaying solutions is to start from an asymptotic distance, 
where that behavior is built-in, and integrate the ODEs inward toward the origin. At large 
distances the ODEs simplify, for this potential, to 

g\x) - xf(x) = , 

f(x)-xg(x) = 0, (15) 
which have solutions at large x = x max , 



/(W) = -a ie - x ^ 2 . (16) 

Note that this asymptotic behavior is the same as that of the non-relativistic simple harmonic 
oscillator wave functions. [3] These forms will be used as starting values (BCs) for the inward 
integration of the full ODEs from x max . Again, ai can be taken as arbitrary for now, to 



be later fixed by the normalization condition, but we need to make an initial guess for the 
energy eigenvalue E. The integration will go from x max back to x — e, again to avoid the 
singularity at x = 0. 

The problem with this shooting inwards method is much the same as that for shooting 
outwards. Unless one somehow is able to guess the exact value of E, the solutions will blow 
up exponentially as one nears the origin, and usually well before that. 

C. Shooting In and Out and Matching 

Thus, after experiencing these well-known problems with ODEs with eigenvalues, we con- 
cluded that we had to combine the two methods, invoking a well-known numerical method, 
"shoot and match." In this case we shoot both outwards from x = e to a point in the middle, 
^match, and inwards from x max back to a; matc h- The two values for g ut(£match) and fi' in (a; m atch) 
will generally differ for a given choice of E and another parameter, which we choose to be 
a . That is likewise the case for /outmatch) and /m(^match)- What we therefore need is 
to find a way to iteratively vary E and ao to reduce the two gaps, Ag and Af, to zero. 
[The asymptotic parameter a\ is chosen to be fixed, i.e., not to be varied, but of a size that 
makes the initial determination of the gaps at x matc h reasonably small; it will in the end be 
determined by the normalization condition, Eq. (fT2|) .] 

We need to define the gaps Ag, etc., at the match point. Initially, we simply chose them 
in terms of their actual differences, 



and likewise for Af, Ag', and Af. Indeed, this often works for the problem we originally 
set out to do, namely, to calculate wave functions for the (nearly) massless u and d quarks. 
However, as the quark mass fh increases, the system becomes more and more non-relativistic. 
This means that the lower component wave function, f(x), becomes relatively small com- 
pared with g(x). Thus a better definition for the gaps would be to make them relative to 
their average value, e.g., 



Ag(E,a Q ) = g nt(x match ) - g in (x mat ch) , 



(17) 



Ag(E,a ) 




(18) 



and likewise for Af(E,ao). 
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D. Zeroing the Gaps 



The two gaps can be driven to zero by using a generalized Newton-Raphson method. 
First, recall how one derives the Newton-Raphson method for finding the zero of a function 
y(x). Make an initial guess for the solution, call it X\. To find a better guess, x 2 , expand in 
a Taylor's expansion, 

y(x2) = y(x 1 ) + y'(x 1 )(x 2 -x 1 )^ — • (19) 

As we want the left-hand side to vanish, the next (and better) guess for the root is 



x 2 = xi- y(xi)/y'(x 1 ) . 



(20) 



This procedure can be repeated until y(x n ) is small enough to be considered zero. 

The generalization of Eqs. ( TT9|) and ( 120]) for our problem is to solve a two-by-two linear 
system for new values of the parameters, E new and a 0i new, from these equations: 

'dAq\ ,~ ~ , fdAg s 



,ncw — «0,oldj 
(oo ,now ^0,oldj 



(21) 



We will calculate the partial derivatives needed above numerically. The procedure to 
find better values of the two parameters E aew and a 0i new is to be iterated until the gaps are 
sufficiently small. 

To solve this linear system in Eq. f[2"Tj) . it is convenient to define a matrix 



M 



d Ag 8 Ag 

dE dao 

dAf d Af 

dE da 



(22) 



and a column vector containing the gaps and partials obtained using the old parameters 

E \,i an d CJo,oldj 

' ( ^id + (^) a old - A, old 
(^)^id+(^f) a old -A /old 
Multiplying C i d by the inverse matrix M _1 then provides us with a column vector containing 
the improved parameters for the next iteration: 



a 



old 



(23) 



M Coid — [E new , ^O.newl 



(24) 
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It turns out that this iterative procedure works very well, with the gaps in the slopes of 
g(x) and f(x) at the match point also going to zero. 

At this point we should remark that others might have done this calculation but by zeroing 
the gaps in the logarithmic derivatives, such as(dg(x) / dx) / g(x) . This has the advantage of 
removing the scale dependence of the functions g and /, but our choice of defining the gaps 
as relative, as in Eq. (JT8]) . is essentially equivalent, as it also involves (semi-)local differences 
divided by the (semi-)local value. 



E. The Is Radial Wave Functions and Some Programming Details 

We implemented the iterative numerical process discussed above by developing Mathe- 
matica notebooks, (l9j but any reasonable programming language could be used instead. 
An example notebook is available from our group's web site. [20] (We especially invite our 
Canadian colleagues to convert this to a Maple program!) Here, as we describe the calcu- 
lation of the ground state radial wave functions for the GMSS potential of Eq. we will 
comment on a number of the programming issues that we encountered. 

The first thing to be done is to set the value of k in the equations to be solved. It is an 
integer and depends upon the angular momentum quantum numbers, Eq. fl5]). For the Is, 
j = 1/2 ground state, I — 0, V — 1, and thus k = —1. 

As noted earlier, we work in units where h = c = 1 {21I so that the dimensions of Eqs. ([7]) 
and ([8]) are fm . For the cases discussed here and below, we have fixed the constants in the 
GMSS potential at k 2 = 0.9 GeV/fm (i.e., « = 2.14 fin" 1 ), and r = 0.705 fm. Making the 
equations dimensionless, as mentioned after Eq. ( ITT]) , this linear scalar potential simplifies 
to V s (x) = x — xq, where x = nr = 1.506. 

For dimensionless distances we chose e = 10~ 6 , x matc h = 1-0, and x max = 6.0. We also 
define a small increment 5 = 0.0001, which will be used when we calculate the partial 
derivatives in Eq. (T22]) numerically and for testing when the gaps to be zeroed are small 
enough. 

The two integrations, inwards and outwards, were done using Mathematica's NDSolve 



function, but one could use any standard Runge-Kutta procedure. [181] To proceed, we 



need the four boundary conditions (BCs) from Eq. ffTS]) and Eq. (fT6"|) . We have defined 
a subroutine, shootinandout, to do this. This subroutine is, of course, a function of the 
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parameters E and ao- It saves the results of the integrations as ins o In and out so In, 
respectively. 

In this case, we already had a good idea Q of what the energy eigenvalue E is for the Is 
ground state. After some fiddling, we found an initial choice of parameters 

E = 0.82, a = 0.2, ai = 1000.0 (25) 

for which the resulting outwards and inwards integrations yielded the curves for g(x) and 
f(x) shown in Fig. [TJ (As mentioned above, the asymptotic normalization a\ is at this point 
arbitrary, to be fixed later by the normalization condition.) 

It is useful at this point to define a function, calcmatchgaps, that uses insoln and 
outsoln to calculate and print out the values of g(x) and f(x) at the match point x matc h, 
along with their slopes and their gaps scaled as in Eq. ( {TBI) . These numbers provide guidance 
as to how to proceed when developing the Mathematica notebook. 

Next, we need the four partial derivatives in the two-by-two matrix M of Eq. (|22"1) . In 
pseudocode, the subroutine for calculating d Ag/dE is 

subroutine dggapbydE(tildeE,aO) 
delE = delta*tildeE 

outsolnl = outward integration with parameters (tildeE + delE.aO) 

from epsilon to xmatch 
g_outl = g(x) at x_match from outsolnl 

insolnl = inward integration with parameters (tildeE + delE,aO) 

from xmax to xmatch 
g_inl = g(x) at x_match from insolnl 
ggapl = 2.0*(g_outl - g_inl)/(g_outl + g_inl) 
dggapdE = (ggapl - ggap)/delE 
return dggapdE 
end subroutine 

Similar subroutines are also implemented for the three other partial derivatives. We 
also found it useful, once these subroutines were in place, to define another subroutine, 
gapsandpartials, which calls shootinandout for the present values of the parameters E 
and ao, followed by calls of calcmatchgaps and the four subroutines for the partial deriva- 
tives. 

At this point we are ready to solve the two-by-two linear system for improved values of 
the parameters E and a . Manipulating arrays is a bit tricky in a programming language 
such as Fortran or C, but in Mathematica one can simply define the matrix M in Eq. ( )22l) 
as a list of lists and the column vector C id is a simple list. Mathematica also provides a 
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built-in function to invert the matrix M, so the evaluation of the right-hand-side of Eq. ( I24j) 
is easy. This procedure is encapsulated in a subroutine we called solvelinsys () . 

Before going on to the program for the iterative loop we need to decide when the gaps are 
small enough to stop the iteration. Thus we defined a test function which returns a boolean 
True if any of the four gaps, Ag and Af as well as their slopes Ag', and Af, is greater 
than 5. Again, in pseudocode, 

subroutine testgapsO 

boolean pi = Abs [ggap] > delta 
boolean p2 = Abs [f gap] > delta 
boolean p3 = Abs [gpgap] > delta 
boolean p4 = Abs [f pgap] > delta 
boolean q = pi OR p2 OR p3 OR p4 
return q 
end subroutine 

The iteration will stop when q = False. 

The coding for running the iteration from an initial choice of parameters E and ao is as 

follows: 

{newtildeE,newaO} = {start_tildeE, start_aO} 
print starting parameters {newtildeE, newaO, al} 
iter = 
q = True 
do loop 

iter = iter+1 

{oldtildeE, oldaO}- = {newtildeE,newaO}- 
shootinandout (oldtildeE, oldaO) 
gapsandpartials (oldtildeE, oldaO) 

solvelinsys () // returns improved values of newtildeE and newaO 
print iter and gaps 
print newtildeE and newaO 
if q 

continue to next iteration, i.e., go to the top of the do loop 
else 

stop and break out of the do loop 
end if 
end do loop 

For the starting parameters used above, E = 0.82 and ao = 0.2, we find that our program 
closes the gaps in four iterations. The final parameter values are E = 0.727102 and ao = 
0.194709. That E correponds, in more conventional units, to an energy of E = 0.306 GeV. 
Figure 2 shows plots of the final, converged g(x) and f(x) after normalization. Note the 
relatively large size of the lower component f(x), showing the importance of relativity in 
this case of massless quarks. 
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Figure 3 shows the corresponding plots for if) a (x) and —ipb(x). As one ought expect for a 
ground state wave function, the upper component if) a (x) has no nodes and the p-wave lower 
component ipb( x ) has one (at the origin). 

We have plotted the t/>'s in Fig. 3 this way to compare with GMSS's Fig. 2. The ^'s 
here are for fitting the non-strange qq mesons and extend out further than those of GMSS. 
The reason for that is that GMSS calculated the u and d quark wave functions to fit the 
non-strange baryon spectrum, which is why they used a value of r = 0.57. Physically, the 
effective origin of the confining potential for mesons is not as strongly localized as in the 
baryon case, where one more quark damps fluctuations additionally. 

F. Some Excited States 

The procedure outlined in the previous sub-section can be applied to calculate the energy 
eigenvalues and wave functions for excited states. One expects that the first excited state 
is the 2s | state. This state also has I — 0, I' — 1 and k = —1, so the BCs at x = e 
are the same as for the Is case. The upper component wave function g(x) should now 
have a new node, i.e., cross the x-axis somewhere between the origin and infinity. Thus, if 
it is desirable to have g(x) start off from the origin going positive, one should choose the 
asymptotic normalization parameter a\ to be negative. 

In view of the asymptotic behavior here being similar to that of a simple harmonic 
oscillator, we might expect that the energy eigenvalue of this state is roughly twice that of 
the ground state. Choosing the parameters needed for the initial integrations to be 

E = 2.1, a = 0.3, ai = -4000.0 , (26) 

the iterative loop again converges in four iterations. The final parameters are E = 1.91897 
(i.e., 0.809 GeV) and a = 0.378981. Plots of ij) a (x) and -ifo(x) are displayed in Fig. 4. For 
the 2s 1/2 state, ip a (x) has one node and ipb(x) has two. 

Why are we counting nodes, anyway? Because a state with more nodes has more energy. 
Qualitatively, more nodes means more curvature, and that means a bigger contribution from 
V 2 , which in turn means a bigger kinetic energy. 

Besides the 2s \ state, there are also two nearby p-wave excited states. The 2p | state 
also has j — I + \, but now with I — 1, I' — 2, and thus k = —2. The BCs near the origin 
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are given by Eqs. (fl3"|) and ffl4l) . but otherwise the coding is very similar to (basically can 
be copied from) the Is | and 2s | cases. Starting this time with 

E = 2.1, a = 3.0, ai = -4000.0 , (27) 

the gaps close in four iterations, with final parameters E = 2.23003 (i.e., 0.940 GeV) and 
a = 1.57784. The 2p | wave functions i/j a ( x ) and — i>b{ x ) are shown in Fig. 5. In contrast to 
the 2s case, here both ipa{x) and ipb{x) have two nodes, reflecting a higher energy eigenvalue. 

The other p-wave state, 2p ~, is qualitatively different. This is the first case for which 
j — I — 1/2, and k is now a positive integer. In a sense, this switches the roles of g(x) and 
f(x). The boundary conditions for the outward integration when j — I — | are different 
from those of Eq. (TlBl : 

/(e) = a e l , g(e) = b e l+l . (28) 
Substituting g(e) into the first ODE allows us to determine this bo, so 

g{e) =a (E + V s (e) + m)e l+1 /(2l + 1) . (29) 

That is, if > 0, both g(x) and f(x) start out from the origin with positive slopes. 
Starting for this 2p ~ case with initial parameters 

E = 2.1, a = 3.0, ai = -4000.0 , (30) 

the iterative procedure converges on the fifth iteration, with final parameter values E = 
2.37846 (i.e., 1.002 GeV) and a = 0.171483. The ^'s are plotted in Fig. 6. Here also both 
component wave functions have two nodes. 

We note in passing that the higher spin state of this pair has the lower energy, contrary 
to the well-known case for the Coulomb potential. This is a common feature for nuclear 
states (even ground states) and reflects the presence of an effective scalar potential as found 



in many nuclear potential models. 



23] 



The energy difference between the 2p | and 2p | states, 62 MeV, is due to a spin- 
orbit interaction. However, there is evidence in the meson spectrum that the spin-orbit 



interaction is suppressed. 24| Page, Goldman, and Ginocchio (PGG) 25j claim this reflects 
a cancellation between a scalar potential V s and a vector potential, V v , having the same 
linear slope at large distances. We will return to this point in Sec. IIVI 
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G. Solutions for Massive Quarks 



The cases for the linear scalar potential V s (r) of Eq. (]9]) when the quarks are massive are 
computed st raig htforwardly using the program described above. For example, using masses 
appropriate 22] for the Qq mesonic states, where q stands for a massless non-strange quark 
[u and d) and Q for the massive strange (s), charmed (c), and bottom (6) quarks, we find 
the following eigenenergies: 



m u = md = 0.0 GeV , E u = E d = 0.306 GeV , (31) 

m s = 0.3445 GeV , E s = 0.486 GeV , (32) 

m c = 1.803 GeV , E c = 1.667 GeV , (33) 

m b = 5.298 GeV , E b = 5.007 GeV , (34) 

by requiring a match to the appropriately weighted average of the lowest pseudoscalar and 
vector states. 26j Note that for the heaviest quarks the states are conventionally bound, 
Eq < rriQ, and confinement need not be invoked to understand the stability of the state, in 
contrast to the situation for the light quarks. 

Figure [7] compares the results for the u, s, c, and b quarks for the upper components, 
ip a {x), and Fig. [8]does the same for the lower components, —ip^x). Note that as the mass m q 
increases, the Qq system becomes more and more non-relativistic, i.e., the lower component 
wave function ipb gets smaller, relative to the upper component ip a . Also, as m q gets larger, 
the wave functions are more and more concentrated near the origin. 

III. VECTOR POTENTIAL HYDROGEN-LIKE ATOMS 

In contrast with the case of massless quarks, an electron bound in a Coulomb potential, 
as in the hydrogen atom, is at the opposite extreme, i.e., is to a very good approximation 
a completely non-relativistic system. The binding energy of the ground state, Ry = 13.6 
eV, is very small compared to the mass of the electron, m e = 0.511 MeV. (In this section 
we will use MeV instead of GeV.) The Schrodinger equation for this problem predicts its 
energy levels basically correctly, missing only the fine-structure splitting, a spin-orbit effect 



about 10 4 times smaller than the binding of the n = 2 levels. 271 ] . 
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As noted in the Introduction, the analytic solutions of the Dirac equation for the hydrogen 



atom were found long ago, [8] and these solutions do give the fine-structure splitting as a 
relativistic effect. Nonetheless, it is an interesting numerical exercise to see if the program 
for the linear scalar potential outlined above (or something like it) can be applied to the 
hydrogen atom case, and with enough accuracy. We want to solve Eqs. ([3]) and (j3J) with 
V v (r) = —Za/r. We will set Z = 1, as these solutions are somewhat more delicate. Heavier, 
hydrogen-like atoms with Z > 1 can be done in a similar manner. 

Converting to dimensionless equations by dividing through by m e (instead of k as before), 

k 

g'{x) + -g{x) - {1 + E + Za/x) f(x) = , (35) 

fix) - -f{x) - (1 - E - Za/x) g(x) = , (36) 

x 

where now x = m e r and E = E/m e . The energy E is better written as E = m e — B, where 
B is the binding energy of the level of interest (which has principal quantum number n and 
orbital angular momentum quantum number I). Note that, because of the smallness of the 
binding energies, E is only slightly less than m e . Thus E is less than (but very close to) 
1, which is why we factored out a minus sign from the third term in Eq. (136]) . relative to 
that in Eq. 01]). We will use the binding energy B = B/m e instead of E as one of the two 
parameters to be determined in the matching of the inwards and outwards integrations. 

A. Boundary Conditions for Inward Integration 

For asymptotically large x the ODEs reduce to 

g'(x) = (1 + E)f(x) , 

f\x) = {1-E)g{x) . (37) 
The solutions of these asymptotic equations at x — x max are 

ff(x mffl ) = ai e -^ w , 

/(aw) = -oi f e""— , (38) 

which will be used as the BCs (starting values) for the inwards integration with a\ as 
another parameter to be determined later by the normalization condition. The coefficient 
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in the exponential decay is small: 

fi = \ll-E 2 = ^2B-B 2 « \f2~B~ < 0.00730 « a . (39) 

The Coulomb wave functions are quite long-ranged, and are all the more so the smaller the 
binding energy. Equation (138]) also shows that, asymptotically, f(x) is very much smaller 
than g(x). The smallness of the lower component is an indication of how non-relativistic 
the hydrogen atom is. 



B. Boundary Conditions for Outward Integration 

As before, to avoid the singularity at x — 0, we integrate outwards from x — e, where 
e takes an appropriately small value. In this case, one might even consider taking e as 
the radius of the nucleus providing the Coulomb potential (but see Sec. IIIIFI for a better 
approach). For the boundary conditions for integrating outwards from x = e when j '• = 1 + |, 
i.e., I' = I + 1 and k = —(I + 1), we can again assume 

g(e) = a e l+1 , f(e) = b e l+2 . (40) 

Substituting these in the second ODE, we can solve for bo, finding 

f(e) = -«o^ + 0(e 1 ^) . (41) 

Noticeably different is the lowering of the power behavior of f(x) due to the 1/x behavior 
of the potential. Note that, because of the factor of a ~ 1/137, the magnitude of f(x) is 
much smaller than that of g(x), also here near the origin. 

For the j — I — | case, i.e., V — I — 1 and k = I, we let f(x) determine the nature of the 
boundary conditions [as in Eq. ( 128]) ]. i.e., 

/(e) = a e l , ^(e) = b e l+1 . (42) 

Substituting these in the first ODE, we again solve for bo and find 

9(e) = a ^je l + O(e^) . (43) 

Again, for this case, g(x) and f(x) start off with the same slope, as in Eq. (129]) . And again, 
there is a lower power behavior of g(x) due to the 1/x in the potential, as in Eq. (]4"1~|) . 
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C. The Is Ground State 



The natural unit of length for the hydrogen atom problem is the Bohr radius, as = 
Ti 2 /m e a = 0.529 A = 0.529 x 10 5 fm. Thus a natural scale for x is in units of the dimensionless 
Bohr radius, xb = aBm e /Tic which isl/ain^ = c= l units. 

In view of the slow exponential decay in Eq. ( 138]) (i.e., the smallness of //), for our 
calculations we chose x max = 7xb, along with x matc h = 0.5xb- Fixing a\ = 10 (recall, it 
is arbitrary before normalization) and the starting values of the parameters to be tuned, 
B = 12 x 10~ 6 MeV (less than the expected 13.6 eV) and a = 0.015, we found the 2 x 2 
matrix of partial derivatives to be 



M 



-108877 107.348 
-320810 107.494 



(44) 



This shows a big sensitivity of the hydrogen atom wave functions to the choice of binding 
energy. Nonetheless, the iterative process of refining B and proceeds nicely to a solution 
with B = 13.6059 eV and a® = 0.0104306. The wave functions g(x) and f(x) are displayed 
as the solid curves in Figs. 191 and [TUl As expected, f(x) is much smaller than g(x). Their 
shapes differ from those shown for the scalar linear potential shown in Fig. |2j 

Figures M and [TU] also display, as dashed curves, the corresponding ip a (x) and 4>b{x) 
radial wave functions. These are, in fact, both purely decaying exponentials proportional 
to e~^ x . This is not surprising for the upper component, tp a {x), since that is the just what 
the Schrodinger equation predicts for the hydrogen atom ground state. It may not be so 
obvious, however, that the lower component, ipb{%), has the same form. That it has no node 
at the origin (in contrast to Fig. [3} is a consequence of the Za/x potential changing the 
presumed x' +2 -dependence of f(x) near the origin to an x' +1 -dependence. 



D. Some Excited States 



The n = 2 excited states of hydrogen are done in a very similar manner to the Is ground 
state calculation. The major difference is that, from the Balmer formula [271 ]. we expect the 
binding energy B for these states to be about 1/n 2 = 1/4 times the ground state binding 
energy. Also, based on Eq. (I39l) . we expect the wave functions will extend outward about 
twice as far, so, for our calculations we chose x max = 16 a B instead of 7 a B . 
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For the 2s \ state, starting with initial parameters B = 3.2 eV (presumably low), ao = 
0.0007, and (fixed) a\ = —10.0, the program converges in four iterations to give final B = 
3.40085 eV and a = 0.00071905. The 2s wave functions are plotted in Figs. HH and Q21 As 
expected on general grounds, unlike the ground state, here i/j a (x) and ipb(x) each have one 
node. 

For the 2p 1 state, starting with initial parameters B = 3.2 eV (also presumably low), 
a = 0.0004, and (fixed) = 100.0, the program converged in four iterations to give final 
B = 3.40086 eV and a = 0.000393546. In fact, the 2s | and 2p | levels are degenerate 



271 ] as a result of an 0(4) symmetry hiding in the equations. The small difference in the 



converged B's we find here is within the numerical precision of our Mathematica program. 
However, the 2s | and 2p | wave functions are very different. The 2p | wave functions are 
shown in Figs. [T3]and[TU Note that, here also, ip a and ipb each have one node. 

For the 2p | state, starting this time with initial parameters B = 3.2 eV, a = 0.000015, 
and ai = 100.0, the program converged in four iterations to give final B = 3.40144 eV and 
ao = 0.0000141876. The difference we find between the j — § and j — \ energy levels, 
0.58 meV, is the spin-orbit splitting, an intrinsically relativistic effect. The analytic value 
of this splitting is 0.453 meV s), [27J and the difference here is due to the limited machine 
precision used in our Mathematica program. The 2p | wave functions are shown in Figs. [15] 
and[16l The ^ a (a;) , s for the 2p | and 2p | states are very similar, but the lower component 
ipb{xys are quite different, reflecting the fact the lower components and the difference are 
both entirely due to relativity. For the 2p | state also, ip a and ipb each have one node. 



E. Hydrogen- like Atoms with Z > 1 

There are no big surprises here, as Table I shows. From Eq. (139 p . the exponential fall-off 
of the wave functions is faster, as ~ \f~B ~ Z, so we adjust the values of x matc h and x max 
accordingly by dividing the hydrogenic values by Z. Running the code for the Is ground 
state for V v (r) = Za/r works well for Z < 100. As Z grows larger, the relativistic corrections 
to the energy become increasingly important. 

As Z approaches 1/a = 137.036, however, the numerical solutions become inaccurate and, 
eventually, unstable. This can be seen in the rapid increase in the value of the parameter 
ao- The value of Z = 1/a is where the Klein paradox [28] comes into play, since the analytic 



18 



result for this state is [12] 



E(ls) = m(l - Z 2 a 2 ) 1/2 (45) 

and the eigenenergy becomes complex beyond that point. The resolution of this paradox 
is, as has been well-known for a long time, the creation of electron-positron pairs from the 



Dirac- Fermi sea when Z > 1/a. 29] 



F. Muonic Atoms 

When a negatively-charged muon slows down in matter, sometimes it is captured by the 
Coulomb potential of an atomic nucleus before it decays. If so, it then quickly cascades 
down through its hydrogen-like levels to the Is ground state, emitting x-rays along the way. 
From there it then either decays or can be be captured by the nucleus, both of which are 
weak interaction processes. The capture process depends sensitively upon the value of the 
Is wave function at the origin. 

The numerical coding for muonic atom ground states for values of Z of interest goes pretty 
much as for the electronic atom case for such a Z. The differences come from replacing the 
mass of the electron (0.511 MeV) with that of the muon (105.66 MeV). This must be done for 
the "muonic Bohr radius", ag^ = as(m e /m M ) = 255.8 fm, and consequently in the choices 
of ^match and x max which are now proportional to as^jZ . Also, the non-relativistic binding 
energy B is now Z 2 (m M /m e ) Ry. 

Calculating, as above, the ground state for muonic 40 Ca [Z = 20), starting with B = 1.1 
MeV and ao = 0.1, the iterations converge to final values of B = 1.13136 MeV and = 



, or a 



0.239462. However, this initial calculation used a value of e = 10 4 (dimensionless' 
cut-off at the origin of 0.000187 fm, rather smaller than the size of the 40 Ca nucleus, 30j 

i? Ca « r A 1/3 = 1.2 x 40 1/3 = 4.1 fm. (46) 

The nuclear charge is not a point charge, of course, but has a charge density distribution, 
p(r), which is spread out over the volume of the nuclear sphere. In fact, one of the major 
reasons for having studied muonic atoms in the past, through the last 2p — > Is x-ray, is for 



determining this charge density, llj 



As a simple example, suppose the nucleus is a sphere of radius R and charge eZ with 
a charge density p which is constant out to its surface and zero for r > R: We need to 
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find the potential seen by the negative muon for r < R. If it were sitting exactly at the 
center of the nucleus, it would (by symmetry) feel no force. That is, there is no term in 
the potential linear in r, as the force at the origin is — dV QUC / dr] r=0 = 0. This suggests that 
we can take V QVLC (r) = A + Br 2 . The coefficient B can be fixed from the requirement that 
the slopes dV nnc /dr and dVconi/dr match at r = R. The coefficient A in turn is fixed by 
the requirement of continuity at r = R. Thus a potential reasonably appropriate for such a 
muonic atom is 

[ aZ (r 2 - 3R 2 )/2R 3 for r < R, 
V v (r) = { \ ~ (47) 

I —aZ/r for r > R. 

This potential no longer has the 1/r singularity at the origin, so the BCs for the outgoing 

integration for the Is state are more like those in Eqs. (|T3|) and ()14p for I = 0: 

g(e) = a e, and /(e) = a (l - E) e 2 /3 , (48) 

bearing in mind that we here are making the equations dimensionless by setting m M = 1 
and that the eigenenergy of the bound state is less than m^. Unlike the pure Coulomb case, 
here /(r) will have curvature near r = 0, where it will vanish. 

Running the numerical integrations with potential V v (r) for 40 Ca gives a binding energy 
B = 1.07006, some 60 keV smaller than the pure Coulomb result for this nucleus. Figure [T71 
displays the difference between the -?/> a (a;)'s and ijj^x)^ for the case with this nuclear charge 
distribution and the pure Coulomb case. 

IV. COMBINING THE TWO POTENTIALS 

We now consider the solution of the radial Dirac equations, Eqs. ([7]) and ([8]), when both 
potentials, the Lorentz scalar V s (x) and Lorentz vector V v (x), come into play. We discuss 
two different situations. 

A. When V v (x) Is Coulombic 

As our first case, we take the scalar potential V s (r) to be the same as in Eq. (jH]) and the 
vector potential as V v (x) = —a s /x. Here a s is the (strong) gluon-quark coupling constant 
of quantum chromodynamics (QCD), having a value of 0(1) at low energies and becoming 
small at high energies (asymptotic freedom). 
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The modifications to the program are minor. The two BC's at x max are the same as in 
the pure V s case, Eq. f|T6|) . The BC's at x = e are slightly modified from that of Eq. f fT3|) 
because here it is V v that dominates. For / = 0, 

g(e) = a e, /(e) = -a (V v (e)/3) e 2 = a a s e/3 , (49) 

as in the hydrogen atom case, but with a s generally much larger than a. 

Running the code for a s = 0.4 gives the wave functions ip a and ipb shown in Fig. [18j For 
comparison, the wave functions for case of a s = (i.e., only using the scalar potential) are 
also shown. The eigenenergies for these two cases are 0.306 GeV (a s = 0) and 0.251 GeV 
(a s = 0.4). That is, adding the Coulomb attraction lowers the ground state energy. 

The sharp rise in ip a {x) and rapid falloff in ipb{x) near the origin is another effect of the 
point Coulomb attraction. This behavior is unphysical, as the motions of the quark and 
anti-quark in the qq meson will smear out the point Coulomb potential somewhat. One 
might think that a Coulomb divergence remains in the relative coordinate, but this is a 
non-relativistic prejudice. Relativistic retardation effects and pair creation and annihilation 
both act to ameliorate this divergence, significantly so in this case, as a s ^> a. This will 
give a potential V v (x) that is similar to the potential used to describe the nuclear size effect 
in muonic atoms, Eq. (14Tj) . We will come back to such a calculation presently. 

Before doing that, however, we present the results of running the massless quarks code 
for a variety of choices for a s . The Is ground state energies and the inital slope parameters 
a are shown in the second and third columns of Table II. As the Coulomb attraction 
increases, the eigenenergy falls off fairly rapidly and eventually goes through zero at a value 
near a s = 1. At this point, like the Klein paradox situation discussed in Sec. IIIIEl it is 



possible to produce qq pairs at liberty. In fact, what is called a "quark condensate" forms, 



a modification of the vacuum state. 31] 



To take into account the smearing of the point Coulomb potential by the motions of 
the quarks in the meson, we again use the potential V v (x) given in Eq. f l4"Tj) . but with 
a replaced by a s and R taken as one half of the electromagnetic radius of the pion, i.e., 



R = | x 0.672 = 0.336 fm. 26( We choose this fraction on the basis that the root-mean-square 
electromagnetic radius represents the average separation of the quark and antiquark, but 
since they are oppositely correlated, each should be at about half of this distance from the 
center of mass of the meson. 
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The boundary conditions (BC's) at the origin for this smeared Coulomb potential (for 
angular momentum j = I + |) now become 

g(e) = a e l+1 , /(e) = -a (E - V v (e) - V s (e) - m) e l+2 /(2l + 3) . (50) 

Running the program with a s = 0.4 and the above R gives curves for ip a (x) and —i/ji,(x) 
with similar upturns and falloffs near the origin as in Fig. [18j but less sharply so. One 
should be aware that, while interesting, these differences in the ijj a ,b{%) near the origin make 
only minor differences in expectation values of operators because of the factor of x 2 in the 
integration over d 3 x; see, e.g., the normalization condition, Eq. (fT2"|) . 

For this smeared potential, the eigenenergies and final parameters for various values of 
a s are also tabulated in Table II. As expected from the above remark, the energy expectation 
value is basically unchanged from the point Coulomb value until one reaches the point where 
the quark condensate takes over. However, the initial slope parameters a do differ from 
those of the point Coulomb case, being less pathological as E goes through zero. 



B. V v as the Cornell Potential and the Spin-orbit Force. 

In this subsection we will assume both potentials V s (x) and V v (x) are linearly confining 
and have the same slopes (i.e., have the same string tension k 2 ). We assume the two 
potentials are displaced from each other with V s (x) lying above V v (x), as in Fig. [19j That 
is, we assume 

V s (x) = x + x s , V v (x) = —a s /x + x — x v , (51) 

or, more appropriately, the smeared version of V v (x) as described above in Sec. IIV Al Here 
x s and x v are parameters to be adjusted to get some desired eigenenergy. They give rise to 
a separation between the two asymptotic potentials of Xd = x s + x v . 

The reason for giving the potentials the same asymptotic slopes is to check that the 



spin-orbit interaction disappears as claimed in Ref. 25j. In doing so, we have changed the 
earlier offset of — xq in V s to +x s since, in GMSS (where only a Lorentz scalar potential 
was used), the negative offset — xq represented a reasonable approximation the one-gluon- 
exchange attraction, roughly an average of the two potentials as shown in the figure. 

Now, however, the BC's at large distances are changed from those in Eq. (|T6|) . since, in 
the equation for g(x), Eq. ([7]), the difference of the dimensionless V s (x) and V v (x) -> i has 

22 



its x-dependence cancel. We need to solve the large-x equations 

g \x) -(E + x d + m) f(x) = 0, / '(z) - 2x g{x) = , (52) 

from which 

g"(x) - 2(E + x d + m) x g(x) = . (53) 



This is a minor modification of the differential equation for the Airy functions. [32J Since we 
require decaying solutions at our starting point x max for the inwards integrations, we discard 
the runaway Bi solution for g(x). The BC's at large x = x max are then 



g{xm^) = ^ Ai [2 1/3 {E + x d + rhY'-x maj , 



f(x max ) = a% 



1/3 



^'[2 1 / 3 (E + x d + m) 1/3 x max ] . (54) 



_{E + Xd + m) 2 _ 

Despite appearances, g(x max ) and f(x max ) still start off with the opposite slopes, as in Eq. 
(USD, since, for t > 0, Ai(t) > and Ai '(t) < 0. 

For calculations, we simplified the displacement between the two potentials by setting 
x v = 0, giving us one less parameter to worry about. The first wave function we want to 
compute is the ground state for the massless quark case, m q = 0, starting out with a s = 
and then adjusting both Xd = x s and a s until we get the desired eigenenergy, noted above 
in Eq. (j3~Ij). 

It turns out that to do so is trickier than in the Vg-only case, presumably because of 
the cancellation between V s {x) and V v (x) in the differential equation for g '(x). Our first 
iterative search with a s = found a highly excited state with E = 7.975, not the desired 
0.7263, for which i/j a {x) had six nodes! After some hand- searching on x s , we found we could 
get a ground state solution, i.e., a ip a (x) with no nodes, looking like the 6-quark ip a (x) of 
Fig. [71 This solution started with x s = 4.0 (very large!) and converged, in six iterations, to 
E = 5.735 (about eight times too large). 

From this a s = solution we were able to march down, in small steps in x s (still with 
a s = 0), to x s = 0.20, with the energy now smaller, E = 2.779, but still a factor of four too 
large. (For smaller values of x s the iterative procedure would not converge, with the matrix 
of Eq. ( |22i) becoming singular.) The wave function ip a {%) for this solution was broader than 
that for the x s = 4.0 case, and the ipb{%) was relatively larger in comparison with i/j a . 

At this point we began marching the value of a s in small steps up from zero, i.e., we 
began turning on the Coulomb attraction. As expected this attraction lowers the energy, 
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and we were able to achieve the desired massless quark eigenenergy of 0.306 GeV with the 
following values: 

x s = 0.20, a s = 1.103, E = 0.7263 . (55) 

The Long Marches did achieve the factor of eight reduction in the ground state eigenenergy. 
The wave functions for this massless quark solution are shown in Fig. [201 Note that ip a {x) 
here begins to resemble the hydrogen- atom ground state ij} a {x) depicted by the dashed curve 
in Fig. [U and it is much narrower than that shown in Fig. [3] . 

The ground state solution says nothing about the PGG claim that, for V s (x) and V v (x) 
having equal slopes, the spin-orbit interaction vanishes. To check this we need to find, at 
a minimum, the eigenenergies of the 2p | and 2p | excited states. For the same values of 
x s and a s as in Eq. ( )55l) and after some searching for good starting values for the iterative 
procedure, we obtained the three n = 2 wave functions displayed in Figs. EU [221 and [231 The 
wave functions here are narrower and more hydrogen-like than in the V s -only case (compare 
with Figs. HI and [6]). The eigenenergies, in GeV, for these states are 

E 2s = 1.020, E 2p i = 0.674, E 2p z = 1.008 . (56) 

The 2p | state is nearly degenerate with the 2s state, some 200 MeV higher than in Sec. 
iHFl but this de generacy is probably accidental. The real surprise here is the huge spin-orbit 
splitting between the 2p \ and 2p § states, not what we expected from the claim made by 



PGG. 
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This difference, it turns out, is due to the presence of the Coulomb attraction term in 
V v {x). If we set a s = 0, then (again for x s = 0.2) we find the following n = 2 eigenenergies 
(in GeV): 

E 2s = 1.757, E 2p i = 1.523, E 2p z = 1.523 . (57) 

The wave functions in this case are similar to those in Figs. ETJ [221 and [23] but are somewhat 
broader and less peaked near the origin. These results show the expected PGG symmetry 
(zero spin-orbit splitting). Note, incidentally, that the 2s state is higher in energy than the 
2p states, as in the experimental charmonium spectrum. [2i| One does conclude, however, 
that having a significant Coulomb contribution in V v (x) will destroy the PGG symmetry. 
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V. SUMMARY 



In this paper we have discussed, pedagogically, three increasingly more sensitive cases for 
solving the radial Dirac equations numerically: 

1. A Lorentz scalar potential, V s , which is linearly confining for quarks. The case of 
massless quarks is necessarily relativistic. As the mass of the quark increases, the 
wave function solutions become more and more nOn-relativistic. 

2. A Lorentz vector potential, V v , whose time-like component is a Coulomb potential. 
We first treated, numerically, the hydrogen atom, which is basically non-relativistic 
problem. We do, however, obtain the relativistic spin-orbit splitting between the 2p^ 
and 2p| states (at the edge of our machine precision). 

3. The combination of scalar and vector potentials, aiming at a relativistic quark model of 
mesons. Of particular interest is when we do (and do not) get a relativistic spin-orbit 
splitting. 

We thank L. J. Curtis, J. L. Friar, C. J. Fontes, M. D. Scadron, and R. Sharma for 
discussions and comments on this work. 
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Figures 
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FIG. 1: The first-pass Is ground state radial wave functions g(x) (solid curve) and f(x) (dashed 
curve) for massless quarks in the GMSS linear potential, Eq. Q, given the initial guess for param- 
eters E and ao of Eq. ([25]) . 
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FIG. 2: The normalized Is ground state radial wave functions g{x) and f(x) for massless quarks 
in the GMSS potential, after iterations have converged to final parameters E and ao- 

28 






FIG. 3: The normalized Is ground state radial wave functions ip a ( x ) (solid curve) and —ip^x) 
(dashed curve) for massless quarks in the GMSS linear potential for fitting qq mesons. 




FIG. 4: The normalized 2s excited state radial wave functions ift a ( x ) (solid curve) and —ip^x) 
(dashed curve) for massless quarks in the GMSS linear potential for fitting qq mesons. 
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FIG. 5: The normalized 2p | excited state ip a {x) (solid curve) and —ipb{x) (dashed curve) for 
massless quarks in the GMSS potential. 
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FIG. 6: The normalized 2p ^ excited state ipa(%) (solid curve) and —ipb(x) (dashed curve) for 
massless quarks in the GMSS potential. Note the considerable differences between these wave 
functions and those in Figs. 4 or 5. 



30 




FIG. 7: The upper component radial wave functions, ip a (x), for massless u and massive s, c, and 
b quarks. 




FIG. 8: The lower component radial wave functions, —ipb(x), for massless u and massive s, c, and 
b quarks. 
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FIG. 9: The normalized hydrogen atom ground state upper component wave functions g(x) (solid 
curve) and tp a (x) = g(x)/x (dashed curve). 
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FIG. 10: The normalized hydrogen atom ground state lower component wave functions f(x) (solid 
curve) and ipb( x ) = f( x )/ x (dashed curve). Note the smallness of these wave functions compared 
to those for the upper component. 
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FIG. 11: The normalized hydrogen atom 2s state upper component wave function g(x) (solid curve) 
and ip a (%) = g{x)/x (dashed curve). 
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FIG. 12: The normalized hydrogen atom 2s state lower component wave function f{x) (solid curve) 
and il>b(%) = f(x)/x (dashed curve). 



33 




0.0001 



0.00008 



0.00006 



0.00004 



0.00002 



- 








- 


/ 


\ 








\ 






1 


\ 




- 

- 


1 

I 


\ 












1 




/ S ^\ 




_ 1 




/ X \ 




1 




/ x \ 








\ \ 




. 1 




\ \ 




. 1 




\ \ 




_ I 




\ \ 




_ I 




\ \ 








\ 




. 1 

_l 




\ >v 

\ X. 




J 




X \ 
X 




J . 








J / 

L / 














1 , ~ - -f - T_ ., 



500 



1000 



1500 



2000 



FIG. 13: The normalized hydrogen atom 2p | state upper component wave functions 0.002 x g(x) 
(solid curve) and ip a (x) = g{x)/x (dashed curve). 
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FIG. 14: The normalized hydrogen atom 2p h state lower component wave functions 0.002 x f(x) 
(solid curve) and ipb(%) = f(x)/x (dashed curve). 
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FIG. 15: The normalized hydrogen atom 2p | state upper component wave functions 0.002 x g(x) 
(solid curve) and ip a (x) = g{x)/x (dashed curve). 
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FIG. 16: The normalized hydrogen atom 2p | state lower component wave functions 0.002 x f(x) 



(solid curve) and ipb(%) = f(x)/x (dashed curve). 
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FIG. 17: The unnormalized Is muonic atom wave functions for 40 Ca (solid curves) compared with 
those for a pure point Coulomb potential with charge Z = 40 (dashed curves). 
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FIG. 18: The unnormalized Is wave functions for massless quarks, tj) a (upper solid curve) and ^ 
(upper dashed curve), for both potentials V v = —a s /r and V s of Eq. ([9]). The lower curves are for 
a s = 0, and are the same as in Fig. The value of a s used for the upper curves was 0.4. 
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FIG. 19: Example of having the (dimensionless) Lorentz scalar V s {x) and (smeared) vector V v (x) 
potentials with the same string tension (i.e., the same slopes at large x.) 
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FIG. 20: The normalized ground state wave functions ip a { x ) (solid curve) and —ifjb(x) (dashed 
curve) when V s (x) and V v (x) have equal asymptotic slopes. This is the solution for the massless 
quark case (u and d quarks) using the parameters given in Eq. (|55|) . Note that the wave functions 
are very narrow, having decayed away by x = 2. Compare this solution for the massless quark case 
shown in Fig. [3l 
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FIG. 21: The normalized first excited 2s wave functions ip a {x) an d —ipb( x ) f° r the case when V s {x) 
and V v {x) have equal asymptotic slopes. Compare this solution for the massless quark case (u and 
d quarks) shown in Fig. 2J Compare also with the 2s hydrogen-atom wave functions, if} a ,b( x ) hi 
Figs, rm and PT21 (dashed curves), noting that the smeared Coulomb attraction in V v now forces 
ipb(0) to be instead of a finite value. 
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FIG. 22: The normalized first excited 2p | wave functions if) a (x) and -ipb{ x ) for the case when 
V s (x) and V v (x) have equal asymptotic slopes. Compare this solution for the massless quark case 
(n and d quarks) shown in Fig. [6J Compare also with the 2p ^ hydrogen-atom wave functions, 
i^afiix) i n Figs. [T3l and [T4l (dashed curves). 
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FIG. 23: The normalized first excited 2p | wave functions tp a {x) and -ipb{ x ) when V s (x) and V v (x) 
have equal asymptotic slopes. Compare this solution for the massless quark case (u and d quarks) 
shown in Fig. [5] Compare also with the 2p | hydrogen-atom wave functions, ip a .b( x ) m Figs. [15] 
and [16] (dashed curves) . 
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Tables 



z 


^ 2 Ry 


D 


(numerical) 


B (analytic) 


10 


1.360 




1.3624 


1.3624 


100 


136.0 




161.6 


161.6 


136 


251.546 




433.8 


448.3 



TABLE I: Numerical results for the Is ground state binding energies B for hydrogen-like atoms 
with nuclear charge Z. The second column displays the non-relativistic eigenenergy and the third 
the result of the numerical integrations. For comparison, the analytic eigenenergies from Eq. (|45p 
are given in the fourth column. Energies here are in keV. 





Point Coulomb 


Smeared Coulomb 


a s 


E (GeV) 


a 


E (GeV) 


a 


0.0 


0.306 


0.195 


0.306 


0.195 


0.2 


0.251 


0.244 


0.251 


0.237 


0.4 


0.194 


0.343 


0.194 


0.305 


0.6 


0.135 


0.551 


0.135 


0.417 


0.8 


0.072 


1.052 


0.074 


0.611 


1.0 


0.004 


2.653 


0.007 


0.977 


1.2 


-.089 


13.46 


-.156 


1.759 



TABLE II: Numerical results from solving the radial Dirac equastions for both V s {x) and a Coulom- 
bic V v (x), showing the variation of the ground state energy as a function of a s . Also shown is how 
the initial slope parameter ao changes with a s . 
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